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ABSTRACT 






Manned aircraft that are intended for surveillance or 
to complete a bombing mission will very likely be engaged 
by surface to-air-missiles having guidance systems based on 
infrared (IR) technology. The objective of this study was to 
characterize via simulation the amount of "cover" that can 
be obtained by dropping from a pre- launched, unmanned 
tactical air launched decoy (TALD) a sequence of pyrophoric 
materials to create an IR cloud, analogous to the 
interference created by microwave chaff, that would protect 
the manned aircraft from the missile. The performance 
analysis is based on a simple reticle based model in which 
the two-dimensional (2D) image is reduced to either a 
composite signal, created by the aircraft, or a composite 
noise, created by the pyrophoric expandable. The analysis 
leads to a computer simulation model producing time and 
space dependent signal-to-noise ratios. It is demonstrated 
that the simulation model can answer questions such as how 
long the materials need to burn, how much intensity is 
needed, what wavelength range is most effective, which 
pyrophoric packets should be dropped, and how many. A 
visual model of the time dependent IR pyrophoric cloud has 
also been created. 
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I . INTRODUCTION 

A. BACKGROUND 

The most effective anti-aircraft weapon ever created 
to-date is the infrared (IR) guided missile. Since their 
introduction into operational use in the early 1950s, IR 
missiles have far exceeded the expectations their designers. 
The enabling technology for practical IR missiles came from 
WWII research on detector materials that could be used to 
fashion detectors sensitive in the IR bands. These initial 
materials were sensitive to radiation in the near-IR (1-2 
micron wavelength) band. Over the last 20 years, advanced 
designs have taken advantage of better detector materials, 
moving the bandpass from the near-IR into the mid-IR (3-5 
micron band) , where engine plume provides a superior signal 
for all-aspect engagement. The three to five micron band 
also has less atmospheric attenuation and clutter factors. 

The need to protect aircraft from attack by effective, 
easily launched IR homing missiles led to the development of 
pyrotechnic flares. The flare needs to have a jamming-to- 
signal (J/S) ratio of greater than one to one versus the 
protected aircraft's signature to lure a missile away from 
the aircraft. Since the observed surface area of the IR 
signature of the flare is quite small, a pyrotechnic flare 
has to operate at very high temperatures to match or exceed 
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the in-band aircraft signature. But high temperatures shift 
the wavelength, and we are forced by the laws of physics to 
increase the size of the flare to provide enough energy at 
the correct wavelength. The flare's intensity is presumed to 
rapidly rise to several orders of magnitude higher than that 
of the target signature and ignite near the target aircraft, 
inside the field of view (FOV) of missile's seeker. 

B . APPROACH 

In this thesis the amount of protection that can be 
obtained by dropping a sequence of pyrophoric materials to 
create an IR chaff cloud will be simulated. A simple 
scenario will be defined in which the packets are dropped 
from a pre-launched unmanned tactical air launched decoy 
(TALD) , and the manned aircraft follows, at a later time, 
close to the same path as the decoy, taking advantage of the 
interfering effect created by the pyrophoric clouds. 

The heat-seeking missile is assumed to be in a 
prelaunch status so that missile aero-dynamics are not part 
of this initial study. The imaging capability of the missile 
is assumed to be defined in terms of a definable number of 
pixels within the known f ield-of-view. Also assumed known is 
the detector sensitivity and detectivity over one or more 
bands of wavelengths. The image of the aircraft is defined 
in terms of a distribution of hot spots. Any interference 
provided by the IR signature of the TALD aircraft are 
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ignored. A visual model for expenditure of the IR chaff 
cloud is created. The sequence and timing for the dropping 
of the packets are assumed to be known. Also, both the speed 
of the manned aircraft and the decoy are assumed known. 

The analysis of the problem can answer questions such 
as how long the materials need to burn, how much intensity 
is needed, what wavelength range and distribution are most 
effective, the rate at which pyrophoric packets should be 
dropped and how many, and lastly, quantitatively what kind 
of signal-to-noise ratio (S/N) deterioration can be 
extracted. The mathematical software MATLAB will be used to 
simulate the pyrophoric material, provide image snapshots in 
time of the pyrophoric material, and calculate the resulting 
S/N. 

Some basic IR theory such as Plank's Law, Wien's Law, 
and the Stef an-Boltzmann Law, are reviewed in Chapter II. 
Radiometric definitions, atmospheric absorption, and general 
characteristics of the different types of IR detectors are 
also discussed in Chapter II. The design of the pyrophoric 
model is discussed in Chapter III. The design of the 
aircraft plume model is addressed in Chapter IV. The 
integration of the pyrophoric and plume model is discussed 
in Chapter V. The first part of the simulation code and some 
S/N curves are presented in Chapter VI. The second part of 
the simulation code and the pyrophoric and plume images are 
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presented in Chapter VII. The use of multiple pyrophoric 
expendables is discussed in Chapter VIII. Finally, 
conclusions from the study of the pyrophoric model and 
future enhancements in modeling are presented in Chapter IX. 
Also, in Appendix A, B and C is a brief description of the 
symbols and the quantities that are used during the analysis 
or in the simulation code, the MATLAB codes, and the 
mathematical derivations, respectively. 
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II . BACKGROUND 



A . INTRODUCTION 

In the year 1800, Sir William Herschel, the royal 
astronomer to the King of England, was conducting an 
experiment with a prism in sunlight. The prism spread the 
sun's rays into a spectrum from violet to red. Herschel 
placed a thermometer in the violet color and recorded the 
temperature. He moved the thermometer through the colors 
from blue to red and noticed that the temperature increased 
progressively. He then moved the thermometer beyond the red 
end of the visible region and the temperature continued to 
increase. Thus, he found energy beyond the red; this energy 
has come to be known as infrared (IR) . 

A portion of the electromagnetic spectrum is indicated 
in Figure 2.1. All electromagnetic radiation obeys similar 
laws of reflection, refraction, diffraction, and 
polarization. The velocity of propagation is the same for 
all. They differ from one another only in wavelength and 
frequency. The portion of the spectrum that includes the 
infrared is depicted in greater detail in the lower part of 
Figure 2.1, where the infrared region is bounded on the 
short-wavelength side by visible light and on the long- 
wavelength side by microwaves. It is convenient to subdivide 
the infrared region into several parts. There are no exact 
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designations for the separation of infrared bands. The most 
universally accepted designations today are as follows: 
visible light up to 3 um is designated near infrared (NIR) , 
middle infrared (MIR) is from 3 to 6 um, and far infrared 
(FIR) is from 6 to 15 um. Beyond 15 um is the extreme 
infrared (XIR) . 



s’ Visit** 



Gamma { 
rays , 

__ J 


; n — 

1 1 

X-rays Ultraviolet 1 * 

i i I 


Infra rad 


! 1 

Radio 

J/EHF SHF UHF VHF " HF 


MF 


IF 


~ vlf\ 


i r 


t ~ i: 


1 .... 1 


I ,1 l 


L □ 




1 1 



aiA 1A 10A 100 A 0.1** 1m 10m 100**^ 0.1cm 1 cm 10cm 1 m 10m 100m 1 km 10 km 100km Wavelength 

3xl0 w 3xl0 1# 3xl0 l7 3xl0 l4 3xl0 l5 3xl0 u 3xl0 u 3xl0 u> 3xl0* 1 3xl0 10 3xl0* 3x10* 3xl0 7 3xl0 4 3x10 s 3xl0 4 3x10? Frequency, Hz 

/ 



n ! 

/ 1 Visible l 


Near infrared 


i i 

j Middle infrared j 


Far infrared 


"T ^ 

J Extreme infrared / 


) j V B G Y 0 R } 




1 




! 


1 1 1 1 1 1 1 1 1 1 1 1 I 1 1 1 1 1 1 1 M 


0.4 0.6 0.8 


1 1.5 2 


3 4 6 


8 10 


15 20 30 Wavelength, u 



Figure 2.1 The Electromagnetic Spectrum. (From Ref [1]) 

B . RADIOMETRIC SYMBOLS 

All objects with nonzero temperature will emit 
electromagnetic radiation. Some objects prove to be better 
radiators than others. The (hypothetical) best radiator is 
one that obeys a set of classical equations called the 
blackbody radiation equations. Some of terms used to 
describe IR sources and radiation are defined in Table 2.1. 
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Table 2.1 Symbols, Descriptions, and Units 



Symbol 


Term 


Description 


Unit 


U 


Radiant 

energy 


Energy transferred by 

electromagnetic waves 


Joule 


P 


Radiant flux 


Rate of transfer of 

radiant energy 


Watt (W) 


W 


Radiant 

emittance 


Radiant flux emitted per 
unit area of a source 


W-cm~ 2 


J 


Radiant 

intensity 


Radiant flux per unit 

solid angle 


W-cm 2 


N 


Radiance 


Radiant flux per unit 

solid angle per unit area 


W -cm ' 2 


H 


Irradiance 


Radiant flux incident per 
unit area 


W ■ cm ' 2 


Wx 


Spectral 

radiant 

emittance 


Radiant emittance per unit 
wavelength interval at a 
particular wavelength 


W ■ cm ~ 2 


8 

i 


Emissivity 


Ratio of radiant emittance 
of a source to that of a 
blackbody at the same 
temperature 


Dimensionless 

J 
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C . BLACKBODY RADIATORS 

The spectral density distribution of the power radiated 
into a hemisphere from a blackbody radiator per unit area of 
the source is given by Planck's Law [2] : 



where 

- Wa. is the spectral radiant emittance of the source with 

units of W-cm'^fxm' 1 

- h is Planck's constant (h=6 . 625xl0‘ 31 joules-s) 

- c is the vacuum velocity of light ( c=3xl0 8 m-s" 1 ) 

- X is the wavelength of the radiation (|im or run) 

- k is Boltzman's constant (k=l . 38xl0~ 23 joules-K’ 1 ) 

- T is the temperature of the source in K 

- c x = 27lhc 2 =3 . 74x10“ W-cm' 2 -|tm“ 

- c 2 = hc/k=l . 439x10“ jxm-K 

Plots of the spectral radiant emittance for a typical range 
of temperatures encountered in our environment are shown in 
Figure 2.2. The spectral radiant emittance is the spectral 
power distribution normalized by the area of the source. 
Note that most of the energy is in the infrared portion of 
the spectrum. As the temperature of the source increases, we 
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note that the area under the curve increases, and the 
location of the peak shifts toward the shorter wavelengths. 




Figure 2.2 Spectral Radiant Emittance vs. Wavelength for 
Blackbodies of Various Temperatures. (From Ref [2]) 



The location of the peak of the spectral radiant emittance 
is found by taking the derivative of Planck's Law with 
respect to wavelength. It can be shown that the result gives 
Wien's Law, the wavelength of the peak spectral radiant 
emittance : 



X. 



2.896X10 3 

T[K) 



[jjm] 



( 11 - 2 ) 
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Knowing the temperature of the source, we can calculate the 
location of the peak emittance and, as shown in Figure 2.2, 
higher temperatures imply lower X a . 

The value of the peak spectral radiant emittance is found by 
substituting (II-2) into (II-l) to obtain 



21.2(2 nhc 2 )T s _ 5 

/ A« -bT 



(II- 3) 



where 



b = 1.286 xlO” I5 W - cm' 1 ■ jUtn' 1 fT\ 



The total emittance W is the total power emitted into a 
hemisphere at all wavelengths per unit area of the source. 
It is found by integrating Planck's Law over all 
wavelengths , 

W = °\w x dA (//— 4 ) 

o 

The result of this integration is 



W = 



2n Vr 4 



I5c z h 



2 7.3 



= oT 



(II- 5) 



where 

G = 5 . 67 x 10 " i2 W ■cm~ 2 °K~ A 

which represents the Stefan-Boltzmann Law and shows that the 
total power per unit area of the source increases as the 
fourth power of the temperature. 
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We frequently pass electromagnetic radiation through a 
spectral filter (like the atmosphere) with a spectral 
response of the form S (X) . The emittance that is passed by 
such a filter is 

W(^,A 2 )= \s(fy x <a (^-6) 

A 

Figure 2.3 is an illustration of one such graph where the 
vertical axis is used to compute the fraction of the total 
emittance that lies between a wavelength of 0 pm and some 
upper wavelength divided by the total emittance over all 
wavelengths . 




Figure 2 . 3 



W(0,X) Relative to Total Emittance. (Ref [2]) 
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D. RADIATION FROM OBJECTS 

In order to calculate the radiation emitted from a 
source, we model the source as either a uniformly radiating 
spherical source (e.g, the sun or a light bulb) or as a 
planar source (e.g, the exhaust port of a jet engine) . 

1. Uniform Spherical Source 

For a spherical source with a known temperature ( in K) , 
an emitting area A, and an emissivity (usually assumed to be 
less than one;i.e., a graybody object [2]), we can 
calculate the radiant emittance of the source as 

w = e&r 4 {ii -i) 



Then we can find the total power from the source as 



P = W-A (//- 8 ) 

and the radiant intensity by realizing that it is emitting 

uniformly into a sphere of 471 steradians : 



J 



P_ 

An 



(II -9) 



2 . Planar Source 

For a planar source, we cannot assume that the 
radiation is symmetrical so the power is radiated 
nonuniformly from such a source. A standard procedure is to 
model the source as a Lambertian emitter, where both power 
and radiant intensity will vary with the observation angle. 
For such a source it can be shown that the radiance is [2] 

_ W _ EOT 4 (// — 10) 

n n 
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E. TRANSMISSION OF IR THROUGH THE EARTH'S ATMOSPHERE 

Most infrared systems must view their targets through 
the earth's atmosphere. Before it reaches the infrared 
sensor, the radiant flux from the target is selectively 
absorbed by several of the atmospheric gases, and scattered 
away from the line of sight by small particles suspended in 
the atmosphere. 



When 


the 


particles are 


small compared 


with the 


wavelength 


of 


the radiation. 


the process is 


known as 



Rayleigh scattering and exhibits a X" 4 dependence. For larger 

particles, the scattering is independent of wavelength. 
Scattering by gas molecules in the atmosphere is, therefore, 
negligibly small for wavelengths longer than 2 urn. Smoke and 
light mist particles are also usually small with respect to 
infrared wavelengths, and infrared radiation can, therefore, 
penetrate further through smoke and mists than visible 
radiation. However, rain, fog particles, and aerosols are 
larger and, consequently, scatter infrared and visible 
radiation to a similar degree. 

In the infrared portion of the spectrum, the absorption 
process poses a far more serious problem than does the 
scattering process. The spectral transmittance measured over 
a 6000 ft horizontal path at sea level is shown in Figure 
2.4. The molecule responsible for each absorption band, 
either water vapor, carbon dioxide, or ozone, is shown in 
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the lower part of the figure. Inspection of Figure 2.4 
reveals the presence of atmospheric windows, i.e, regions of 
reduced atmospheric attenuation. 

IR detection systems are designed to operate in these 
windows. Combinations of detectors and spectral bandpass 
filters are selected to define the operating region to 
conform a window to maximize performance and minimize 
background contributions. 
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Figure 2.4 Transmittance of Atmosphere over one Nautical 
Mile (NM) Sea Level Path (Infrared Region). (From Ref [3]) 
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IR DETECTORS 



F. 

The detector is the heart of every IR system because it 
converts scene radiation into some other measurable form; 
this can be an electrical current, or a change in some 
physical property of the detector. 

In general, an infrared detector in a particular set of 
operating conditions is characterized by two performance 
measures: the responsivity R and the specific detectivity 
D* . The responsivity is the gain of the detector expressed 
in volts of output signal per watt of input signal. The 
specific detectivity is the detector output signal- to-noise 
ratio for one watt of input signal, normalized to a unit 
sensitive detector area and a unit electrical bandwidth. The 
different types of IR detectors may be divided into two 
broad classes, namely thermal detectors and photon, or 
quantum, detectors. 

1 . Thermal Detectors 

Because thermal detectors have been used since 
Herschel's discovery of the infrared portion of the 
spectrum, it is appropriate to consider them first. They are 
distinguished as a class by the observation that the heating 
effect of the incident radiation causes a change in some 
physical property of the detector. 

Since most thermal detectors do not require cooling, 
they have found almost universal acceptance in certain field 
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applications in which it is impractical to provide such 
cooling. Because (theoretically) they respond equally to all 
wavelengths, thermal detectors are often used in 
radiometers. The time constant of a thermal detector is 
usually a few milliseconds or longer, so that they are 
rarely used in search systems or in any other application in 
which high data rates are required. 

2 . Photon Detectors 

Most photon detectors have a detectivity that is one or 
two orders of magnitude greater than that of thermal 
detectors. This higher detectivity does not come for free, 
however, since many photon detectors will not function 
unless they are cooled to cryogenic temperatures. Because of 
the direct interaction between the incident photons and the 
electrons of the detector material, the response time of 
photon detectors is very short; most have time constants of 
a few microseconds rather than the few milliseconds typical 
of thermal detectors. 

Finally, the spectral response of photon detectors, 
unlike that of thermal detectors, varies with wavelength. 
Curves of D* versus wavelength are shown in Figure 2.5 for 
different types of photon detectors. 
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Figure 2.5 D* vs. X for Representative Detectors (From Ref 
( 2 ] ) 
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III. PYROPHORIC MODEL 



A. GENERAL DESIGN REQUIREMENTS 

Infrared decoys are used in modern times to protect 
combat systems from IR tracking threats. The most common 
examples of IR decoys are IR flares. IR decoys are used to 
protect aircraft from heat-seeking missiles. Some key 
requirements for the design specifications of an expendable 
decoy are discussed in the following. 

1. Peak Intensity 

Peak intensity is normally the most important 
requirement. IR decoys must radiate with sufficient 
intensity and at least exceed the intended target's radiant 
intensity in the band of interest. This can be accomplished 
by controlling the decoy temperature or by the use of 
selective emitters that have a higher emissivity in the band 
of interest. A common design problem is illustrated in the 
spectral intensity versus wavelength plot of Figure 3.1. The 
relative spectral distribution between a decoy (a small hot 
source) and an aircraft target (a large, relatively cool 
source) is shown. From the viewpoint of the decoy designer, 
it is advantageous if a significantly higher percentage of 
the optical signal at the detector, within the ' decoy band' 
is due to the decoy rather than the target. Finally, the 
peak intensity is the primary driver for the decoy weight, 
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volume, and cost that must be smaller, lighter, and 
cheaper , respectively , than the protected target. 




Figure 3.1 Typical Decoy and Target Spectra (From Ref [4]) 



2 . Flare Intensity Rise Time 

A decoy must persist long enough to ensure that no 
possibility of target reacquisition remains. Consequently, 
it must maintain a credible signature until the original 
target is no longer in the threat f ield-of-view. If this is 
not the case, it is necessary to deploy a second flare. 
Additionally, the decoy must achieve an effective intensity 
quickly enough to capture the seeker before leaving the 
threat f ield-of-view . The diameter of the threat field-of- 
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view at the time of decoy deployment is usually less than 
200 m. This means that effective operational intensity 
levels must be achieved in a fraction of a second [Ref 4]. 
The exact value of the time to rise to peak radiant 
intensity depends of the chemical composition and packing. 
Exact values are difficult to find in the open literature. 

B. DEVELOPMENT OF THE PYROPHORIC MODEL 

The pyrophoric model must simulate changes in both the 
time and space domains. The modeling procedure for 
calculating the in-band radiant intensity is summarized in 
Figure 3.2. 



/ Specify 

Pyrophoric 

Inputs 




Figure 3.2 Calculation of In-Band Radiant Intensity. 
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1. Time Domain Variations 

The emittance over all wavelengths for any hot object 
follows Planks radiation equation is given from the Stefan- 
Boltzmann Law 



where the Stefan Boltzmann radiation constant is [Ref 2] 

<j = 5.67x10 ~ U W ■ cm' 2 ° K~* 

£ is the emissivity of the pyrophoric and T pf a specified 
temperature. The approximation of a graybody model is used 
to approximate the in-band radiant emittance for these 
expendables [Ref 4] . 

The M-file plankpf.m recorded for reference in Appendix 
B is used to calculate for a particular temperature the 
integrated emittance over a band of wavelengths. The 3-to 5- 
micron band is very common for a detector because this band 
has less atmospheric attenuation. The emittance W(X 1/ X 2 ) p£ 
that lies inside the detector's band is obtained using the 
M-file program cited above. The fraction of the emitted 
power that lies between two wavelengths X ; and X 2 is 



The pyrophoric source is modeled as a uniformly 
radiating spherical source with radius r^ at the instant 



W{A,T) P/ = eoT pf 4 



(ill - 1 ) 




(777-2) 
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for which a maximum occurs in the radiant intensity t . The 
associated area is 




(HI- 3 ) 



The peak radiant intensity of an isotropic radiator is 

[Ref 2] 



The total radiant intensity that lies inside the detector's 
band i s 



From the input t p for the pyrophoric burn, the time 
radiant intensity distribution of the pyrophoric burn can be 
predicted. Two operational properties are desirable. 
Ideally, the pyrophoric time function will achieve peak 
intensity quick enough to capture missile seeker attention 
before leaving the f ield-of-view. Furthermore, it should 
persist long enough to ensure that no possibility of target 
reacquisition remains. The probability density function 
(pdf) of the standard Gamma distribution [Ref 5] has 
sufficient flexibility to model these two basic properties. 
The flow chart that summarizes the development of the 




(in-*) 




("I- 5 ) 
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temporal variation in radiant intensity distribution of the 
pyrophoric burn is shown in Figure 3.3. 



Set 

Pyrophoric 

Inputs 



7 




C Create the time radiant intensity 
distribution of pyrophoric bum 




Figure 3.3 The Time Radiant Intensity 
Structure of Pyrophoric Burn. 



Distribution 
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Other inputs such as that corresponding to initial time t t 
and final time t 2 for the simulation to be performed will 
lead to the development of a time-function-array t. 

The Gamma distribution, which is a continuous random 
variable (RV) , is defined as 



t> 0 



P a r{a) 



t ^ e % 



a > 0 , P > 0 



(III -6) 



Because of interest in modeling the temporal variation in 
radiant intensity, the t parameter is used to define time 
dependence. The standard Gamma distribution [Ref 5] has (3=1 
so the pdf of a RV modeled as a standard Gamma RV is 



(™- 7 ) 

r >0 

a > 0 

In Appendix C, a short proof is presented that demonstrates 
that the peak in the distribution satisfies the condition 

1 p = < 3-1 (///- 8 ) 

From this relation the values for the parameter a for 
different values of the pyrophoric peak times are 
predicted. Using the M-file gammaz.m, recorded for reference 
in Appendix B, we calculate the Gamma function T(a) and 
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employ the main program to predict the temporal variation in 
the radiant intensity as 






(III -9) 



A plot of equation (III-9) for the conditions 
T p£ =2000K, and r MAX =1.0m for a pyrophoric time array 
8s are shown in Figure 3.4. 



t p =lsec / 

lasting 



x io 5 




Figure 3.4 The Time Radiant Intensity Distribution of 
Pyrophoric Burn 
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2 . Space Domain 

The flow chart that is followed to develop 
radiant intensity distribution of pyrophoric burn in 
space domain is shown in Figure 3.5. 




Figure 3.5 Normalized Space Distribution Structure 



the 

the 
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It is assumed that the spatial distribution for each 
instant of the pyrophoric burn is a Gaussian distribution. 
Generally, a continuous random variable X is said to have a 
Gaussian distribution if the pdf is 

f(x; ju, a) = -p=L ~ < * < oo (/// - 10) 

v 2-n -a 

where jj. is the mean and a is the standard deviation. One 
Gaussian distribution is shown in Figure 3.6. The Gaussian 
pdf is symmetric about the mean and the standard deviation 
is the distance from the mean to the point at which the 
slope of the curve changes. 




Figure 3.6 Graph of a Gaussian Distribution 



28 



In Appendix C a proof is presented that demonstrates that 
for a Gaussian random variable Z, the area under the curve 



from zero to infinity is given by equation (C-12), 
reproduced here for convenience: 


which is 




(in- 11 ) 



From the inputs t pf r^, and the pyrophoric time function 



array t, the following parameters can be obtained: 




T 

u p — (cons tan? pyrophoric velocity) 

l p 


{ill- 12) 


’'maxi — u p't ( temporal dependent pyrophoric radius) 


(///- 13) 


]T 

r NORM - MAX ' ( normalized pyrophoric radius) 

r MAX 


(/// — 14) 



At each instant of the pyrophoric burn, it is assumed that 



the spatial radiant intensity function is a 
distribution with a mean value 


Gaussian 


_ ^MAxt _ u p't _ t / 
r MAX U p ' tp / P 


{ill- 15) 


From (III-ll), the area under each curve from 
infinity is 


zero to 
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(/// — 16 ) 




which is positive and less than one. 

For each element of the array t, the value of the 
temporal distribution is given by 

/(' ,<! ) = rfr'*"' ' e ~' (m-n) 

I» 

Now in equation (III-10) we use x= r N0RM and (i= t/t p to obtain 

the normalized space distribution in the pyrophoric burn 
as 



f(f, 



NORM ’ 





(hi — is) 



The next step is to use 1/Q(t) to normalize the corresponding 
equation (III-18). The result, which is a pdf defined over 
positive r N0RM , has an area under the curve equal to one. 
This result is multiplied by equation (III-9) to obtain 



J 



_ , \ 

t 



r NORM ’ / 1 



~f( r NORM’(/^ jO’)Xj(t)X( J/ q(J ^ 



(ill- 19 ) 



where the space integration of the space-time distribution 
satisfies the relation 
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This forces the resulting space-time distribution j to have 
a space integrated area equal to J(t) . One example is shown 
in Figure 3.7 where successive curves represent time steps 
of 0.2s. 



x 10 




Figure 3.7 Pyrophoric Space Distribution 
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The last major step, see the flowchart in Figure 3.5, is to 
create the space domain pyrophoric integrated view-model. 
The general idea is to pick a curve from Figure 3.7 which 
corresponds to one instant of time t, as illustrated in 
Figure 3.8. 




Figure 3.8 Pyrophoric Space distribution for one instant of 
time t . 

As represented on Figure 3.9, the pyrophoric source is 
modeled as spherical. It is assumed that an external remote 
observer will 'view' the source and perceive radiation 
integrated along a thin line x. The total radiant intensity 
distribution at a specific instant of time will consist of 
the integration of radiant intensity over each distance | rj, 
where | rj is the distance between the center of the 
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distribution and some point along the line x. The line x is 



shifted by | rj from the distribution center. 




Figure 3.9 The integrated pyrophoric model projection 

Analytically, for each distance from the center of the 
distribution | rj, a fixed number of uniformly spaced points 
are selected along the line x from x=0 to x=x MAX . The 
distance | rj that connects each of these points with the 
center of the distribution is given by 

(//, ~ 20) 

where | x.| the distance along the line x. In the computer 
code, an index k is generated as 
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k = 




\ 




x size{r NORM ) 



(/// — 21 ) 



so the integrated-view radiant intensity for the specific 
| rj, i.e, one point is given by j t (k) . The summation of all 

the radiant intensity values for each | r.| give the 
integrated-view radiant intensity 

J v{r 0 ,t) = Yj X "* X j,( k ) k index (ill -22) 

which is easily computed. Alternatively, an equivalent but 
more analytical representation is 









Jr 






r i’ 




(ill - 23) 



Typical distributions of integrated-view radiant intensities 
for pyrophoric burns lasting eight seconds are shown in 
Figure 3.10. The horizontal axis of Figure 3.10 is specified 
in terms of the index variable of the array r o . It is 
observed that the curve with the maximum area corresponds to 
the instant t=t p . 
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Figure 3.10 Integrated-view Radiant Intensity Distribution 
Curves for the Pyrophoric Model. 
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IV. PLUME MODEL 



A. GENERAL DESIGN REQUIREMENTS 

The plume of the aircraft can dominate the IR 
signature. This is no doubt due to the considerable radiant 
energy created from the combustion process. Measurements of 
the radiation from military a aircraft's plume are highly 
classified. Fortunately, it is not too difficult to make 
order-of-magnitude calculations of the radiation by using 
only temperature and dimensional information that can be 
found in the open literature [Ref 4] . 

During afterburning, which results in an increase in 
thrust, the rate at which fuel is consumed increases 
drastically and the plume becomes the dominant source. Also, 
if the aircraft is viewed from the forward hemisphere or 
from an aspect angle at which the tailpipes are not visible, 
the plume is the only source of radiation available. The 
calculation of the exact radiant intensity of the plume is 
extremely difficult since both temperature and emissivity 
vary in a complex manner through out its volume. 

The exhaust temperature contours for a typical turbojet 
engine with and without afterburner are shown in Figure 4.1. 
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Figure 4.1 Exhaust Temperature Contours for a JT4A Turbojet 
Engine with and without Afterburner (From Ref [1] ) 

We can see that when afterburner is turned on, the 
temperature and size of the plume increase appreciably. If 
we integrate over the entire plume in order to estimate its 
radiant intensity, the radiant intensity will be several 
times that of the hot tailpipe. For in-band engineering 
calculations, a plume can be approximated as a graybody with 
an emissivity of 0.9 [Ref 1] . 

The steps that are followed to calculate the radiant 
intensity of the plume model is shown in Figure 4.2. 
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Specify 
Plume Inputs 




Figure 4.2 Calculation of 



In-Band Radiant Intensity. 
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B. DEVELOPMENT OF THE PLUME MODEL 

The assumed geometric model for the plume, shown in 
Figure 4.3, is an ellipse. This is used to model a three- 
dimensional radiating ellipsoidal source. The required 
parameters are the included major axis a pi , the minor axis 
b pl , and the temperature T P1 that is assumed constant 
throughout the plume. The area of the ellipse is 



Figure 4.3 Plume's Radiating Source Area. 

The Stefan-Boltzmann Law is again used to predict the 
radiant emittance. Here it is applied to the case of the 
plume : 



where £ here is the effective emissivity of the plume. 

Using the M-file plankpl.m for the particular 
temperature with the specific wavelengths of the detector's 
band, we can compute the total emittance W(X 1 ,X 2 ) that lies 




(IV - 1 ) 




W(X,T) n =eoT„‘ 



(IV -2) 
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inside the detector's band. Next, we find the fraction of 



the emittance that lies between wavelengths X 1 and X 2 as 



Vm . 



PI 



W(A,T)„ 



(IV -3) 



The physical model for the plume source is assumed to behave 
as a Lambertian emitter [Ref 2] . Following the standard rule 
for calculating radiance from total emittance [Ref 2], we 
get 



N 



pi ~ 



eo-Tpt 

n 



(IV- 4) 



The total radiant intensity is then given by 



ri ri ri , 

(IV- 

and the plume radiant intensity that lies inside the 
detector's band is 



J (A, , /t- ) w J pj • T] Pl 

(IV -6) 

For later reference, the composite image signal associated 
with the plume target is defined 

S = (IV -7) 
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V. 



INTEGRATION OF PYROPHORIC /PLUME MODELS 



The main program, recorded in Appendix B, consists of 
two parts. The general structure of the main program is 
represented in Figure 5.1. 




Generate S/N Ratio with \ 
variation in parameters I 
over ranee of time J 



/ Generate images at 
one instant 






Figure 5.1 General structure of the main program. 



Once the inputs have been specified, a selection in the 
main program permits execution of either the first part, 
which is the generation of the signal- to-noise ratio (S/N) , 
or the second part, which refers to the generation of the 
pyrophoric flare and plume images. 
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The main steps of the procedure which are followed to 



create S/N are shown in Figure 5.2. 




Figure 5.2 S/N Structure. 
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A. CALCULATION OF S/N 

The first step of the flowchart in Figure 5.2, the 
computation of the time radiant intensity distribution of 
the pyrophoric burn, was covered in Chapter III. One typical 
distribution is represented in Figure 3.4. 

In the second step of the flowchart, the in-band 
radiant intensity calculation of the pyrophoric material for 
each instant of time t is performed. This is found by 
performing a multiplication of the temporal variation in the 
radiant intensity J(t), given by the equation (III-9), with 
the fraction of the emitted power t] p£ that lies inside the 
detector's band (X 1 X 2 ) , given by equation (III-2). 

Next, the in-band radiant intensity of the pyrophoric 
material for each instant of time t is obtained from 

j(A l ,A 2 ) P f-j(t) -Tjpjj (y — i) 



and 



w = y(i,,i 2 ), / (y-2) 

where N is interpreted as a composite noise interference for 
the missile tracking system. 
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The in-band radiant intensity calculation of the plume, 
as given by equation (IV-7) in Chapter IV, is the last 
component needed in the S/N prediction. 

For the calculation of S/N, a situation similar to that 
shown in Figure 5.3 is assumed, where the plume represents 
the signal and the pyrophoric the noise. 



(signal) 



(noise) 

Figure 5.3 S/N Representation. 

Based on the preceding discussion, by taking the ratio of 
the result of equation (IV-7), which represents the 
composite plume image signal, with the result of equation 
(V-2), which represents the composite pyrophoric image 
signal, we obtain an estimate for S/N: 

S _ )pi /y _ 

N " J(\,H 2 ) P/ 
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After the estimation of S/N, the main program permits a 
change in parameters such as pyrophoric temperature T p£ , 
plume temperature T P1 , detector's band wavelengths (X 1 ,X 2 ) , 
pyrophoric peak time t , or radius r^, and a recalculation 
with the new inputs produces a new value for S/N. 

The second part in the main program refers to the 
generation of the pyrophoric and plume images. The main 
steps of the procedure which are required to create the 
images are shown in Figure 5.4. 




Figure 5.4 Generate Images Structure. 
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B. GENERATION OF IMAGES 

1. Pyrophoric Image Construction 

The first four steps of the flowchart in Figure 5.4 
have been already created. Specifically, this includes the 
in-band radiant intensity calculation for each instant of 
time t. This calculation is done in the first part of the 
main program which also deals with the S/N generation. Also 
covered in these four steps is the conversion from the 
normalized space distribution to the integrated pyrophoric 
model discussed in Chapter III. 

To construct the pyrophoric image, it is assumed that a 
situation similar to Figure 5.5 is applicable, where FOV 
represents the f ield-of -view of the incoming missile, R is 
the distance at which the missile will detect the pyrophoric 
source, and N pix the number of pixels that have to be used in 
order to visualize the pyrophoric image. These quantities 
are some of the inputs which are specified at the beginning 
of the main program. 

From this assumption, the missile's window dimension 
Rj. ield is calculated as 

R Field = FOV{rad)x R(m ) [m] (V' - 4) 
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To specify the distance of the pixels from the center of the 
pyrophoric image, the Pythagorean relation is used 






2 + j Z 



(V-S) 



where the i, j indices take values from _ (N pix /2) to (N pix /2). 
Once N r has been calculated for all pixels, the 
corresponding radial distances are predicted from 




N p u 



(V-6) 
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Next, in the algorithm an index is generated in order to 
assign a value from the integrated-view radiant intensity, 
equation (III-23), to each pixel. This index is determined 
by 



For each value of this index that corresponds to the inside 
of the integrated-view radiant intensity array index, the 
corresponding radiant intensity value is assigned 



otherwise the zero value is given to the pixel. Finally, a 
normalization is performed. The value of each pixel is 
divided by the summation of all the pixels and multiplied 
with the in-band pyrophoric radiant intensity, equation (V- 
1) , ( , X 2 ) pe - The normalization rule is therefore equivalent 

to 



i r MAX ■ maX ( r A'ORM )) 



x size(r N0RM ) 



(V-l) 




<y-*) 




(V-9) 



where the time dependence is implicit. 
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2. Plume image construction 

From the flowchart in Figure 5.4, the only necessary 
information to create the basic plume image is the 
calculation of in-band radiant intensity of plume material, 
which has been already estimated. Specifically, the in-band 
radiant intensity calculation for each instant of time t is 
given from equation (IV-7) . 

To construct the plume image, it is assumed that the 
geometric conditions represented by Figure 5.6 are 
applicable . 




Figure 5.6 Plume Image Construction. 
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Note the main difference between Figure 5.5 and Figure 5.6 
is that the plume is represented here as an ellipse instead 
of as a circle. The FOV, R, and N pix are the same quantities 
that already have been specified for the pyrophoric image 
construction. 

To specify the location of the pixels in the plume 
image, one vector r is created 



i,j denote indices that take values from -(N pix /2) to 
(N pix /2) with a specific step 

e x , e y denote unit vectors along the major axis a pl and 
minor axis b. 

Pi 



In order to create an ellipsoidal shape, a value has to be 
assigned to the pixels that belong inside the ellipse's 
boundary 



r =iA-e x + jA ■ e y 



(V-10) 



where 




O' -11) 




(V-12) 



where x,y are defined as 



x = r • e x =iA 



(V-13) 



y = r*e y =jA 
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Taking the preceding into account and making use of 
equations (V-ll) , (V-12), (V-13), we get 




(V-14) 



which provides a convenient rule for defining the interior 
part of the plume in the computer model . 

The numeric values which are assigned to the pixels 

that define the in-band plume radiant intensity J(X 1 ,X 2 }^ 1 are 
taken from plume model and equation (IV- 6) . After the 
generation of the ellipsoidal shape, a normalization takes 
place. The value of each pixel is divided with the summation 
of all the pixels and multiplied with the in-band plume 
radiant intensity J (A, lf X 2 ) P1 . This normalization rule is 
equivalent to 



The similarity of equation (V-15) to equation (V-9) for the 
pyrophoric image is apparent. 




(V-15) 
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VI. PRESENTATION OF S/N CURVES 



As mentioned in Chapter V, the first part of the main 
program, which is recorded in the Appendix B, is primarily- 
concerned with the generation of S/N. Once the inputs have 
been specified, the value of S/N can be extracted for each 
instant of time. A temporal variation in S/N during the 
pyrophoric burn time can then be obtained. Next, by changing 
parameters such as pyrophoric temperature T pf , detector band 
(X 1 ,X 2 ), pyrophoric peak time t p , or radius r^, the program 
computes S/N curves for the new inputs. 

A selection of these curves which have been generated 
through the procedure described are presented in this 
chapter. A graph that presents an overall variation of S/N 
during the pyrophoric burn time and a chart that records the 
minimum S/N for different parameters are included in each 
figure . 

Specifically, Figure 6.1 is a plot of S/N as a function 
of the pyrophoric function time for different values of 
pyrophoric temperature T pf . Results were generated for 
T pl =1000K, t p =lsec, and r MAX =2m. As expected, an increase in 
pyrophoric temperature T p£ results in a decrease in S/N 
during the burn period of the pyrophoric material. This is 
an expected result because it is known from the Stefan- 
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Boltzmann Law that the emittance increases as the fourth 
power of the temperature. 

Figure 6.2 is a plot of S/N as a function of the 

pyrophoric function time for different values of radius r^ 
corresponding to the instant t p for which a maximum occurs 
in the radiant intensity. Results were generated for 
T p£ =2000K, T pl =1000K, and t p =lsec. We observe that increasing 

the radius r,^ significantly decreases S/N between the 
initial and final time of the pyrophoric function. This 
leads to the conclusion that the radiant intensity output of 
the pyrophoric flare increases with r^. 

Figure 6.3 is a plot of S/N as a function of the 

pyrophoric burn time for different detector bands. Results 
were generated for T pf =2000K, T pl =1000K, t p =0.8sec, and r I1AX =2m. 
The curve with the lowest values corresponds to the 1-2 

micron band. This band exhibits a large amount of 

atmospheric attenuation and is not a practical window of 
detection for a missile. The next band is the 3-5 micron 
band, which has less atmospheric attenuation. This band, 
along with the 8-12 micron band, exhibits the lowest amount 
of atmospheric attenuation and are practical windows of 
detection for a missile. Making a comparison between these 
two bands, we see from Figure 6.3 that the curve with the 
lowest values corresponds to the 3-5 micron band. The next 
lowest band is the 7-9 micron band. Its curve has lower 
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values of S/N compared with the 8-12 micron band, but again 
this band exhibits a large amount of atmospheric attenuation 
and is not a practical window of detection for a missile. 

Figure 6.4 is a plot of S/N as a function of the 
pyrophoric function time for different values of fc p . Results 
were generated for T p£ =2000K, T pl =1000K, and r MAX =2m. For this 
case, the minimum S/N corresponding to each t are the same. 
The focus here is the length of time for which the 
pyrophoric flare will achieve intensity that can capture the 
missile's seeker. From this plot, we conclude that the flare 
achieving a peak radiant intensity at a later point in time 
will have an extended windows of effectiveness. This point 
is demonstrated in Figure 6.4 by comparing the length of 
time S/N is below 0.5. The pyrophoric flare reaching its 
peak at 1.4 seconds has a window of effectiveness of 3.7 
seconds, while the flare with t p =0.8 is only effective for 
2.9 seconds. The former is "better" by approximately 25%. 
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Min S/N 




Figure 6.1 S/N for various T pf 
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Min S/N 
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Figure 6.2 S/N for various r^ 
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Figure 6.3 S/N for various Detector Bandpass 
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Figure 6.4 S/N for various t p 
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VII. PRESENTATION OF IMAGES 



The second part of the main program, which is recorded 
in Appendix B, consists of the generation of the pyrophoric 
and plume images. The procedure which is followed is 
described in Chapter V. 

The image visualization is created by using a definable 
number of pixels. The command "image" in MATLAB creates an 
image graphics object by interpreting each element of the 
pyrophoric or the plume matrix which has been created as an 
index into the figure's colormap. Each element of the above 
matrix specifies the color of a rectilinear patch in the 
image. Along with the command "image", the command 
"CdataMapping/scaled" is used to scale the values. 

The variation in radiant intensity for the pyrophoric 
image can be described better with two colormaps. The "jet" 
colormap, which ranges from blue to red and passes through 
the colors cyan, yellow, and orange. The "gray" colormap 
returns a linear grayscale colormap. In Figures 7.2 to 7.6, 
a colorbar is used to show the current color scale. 

The second part of the M-file pyrof.m, which is 
recorded in Appendix B, is used to provide image snapshots 
for different instances of time for the pyrophoric model and 
for the plume model which is constant with time. For example 
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see Figure 7.1 which was generated with T pl =1000K. A sequence 
of pyrophoric images on a "jet" colormap for a pyrophoric 
with a burn function of seven seconds, burn temperature of 
T p£ =2000K, detector band 3-5 micron band, peak radiant 
intensity at t p =1.5 seconds, and radius at t p equal to two 
meters are shown in Figures 7.2 through 7.6 for times t=l, 
t=1.5, t=3, t=5, and t=7 seconds, respectively. Looking at 
these plots, we observe that the ring with the maximum 
radiant intensity will be small and close to the center for 
small times. At the end of the burn, for the pyrophoric 
flare the circular periphery of peak radiant intensity has a 
larger size but a lower radiant intensity value. 

One byproduct of the time-space modeling of the 
pyrophoric burn is that the radiant energy profile does not 
remain peaked at the center. This point was previously made 
in Chapter III and qualitatively represented with the time 
dependent graph of Figure 3.8. The images of the pyrophoric 
flare seen in Figures 7.2 to 7.6 demonstrate this property. 
This is a byproduct of the assumption that the peak in the 
radiant intensity will occur for t p >0 for a spherical "shell" 
at r=r MAX >0 . This is consistent with a physical model of an 
expanding material for which the center burns first and then 
an outward radial wave of combustion follows. 
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Figure 7.1 The Plume Image. 
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Figure 7.2 The pyrophoric image at t=l sec. 
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The pyrophoric image at t=t p =1.5 sec. 
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.4 The pyrophoric image at t=3 sec. 
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Figure 7.5 The pyrophoric image at t=5 sec. 
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Figure 7 . 6 The pyrophoric image at t=7 sec . 
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Representation of the images may give the impression that 
setting the plume in the center of the pyrophoric image will 
allow the missile to detect the plume, which is uncovered as 
the radial wave of combustion moves outward. This is not 
true since the reticle based detection system, represented 
in Figure 7.7, converts all the pixels inside the FOV to one 
value . 




S(,)+N(l) 



Figure 7.7 Reticle based Detection System. 



It should be noted that the S/N predictions in Chapter 
VI are based on this model. In particular, using a reticle, 
the missile electro-optic (EO) process converts two- 
dimensional (2D) optical image to a one-dimensional (ID) 
electronic signal [Ref 1] . 
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VIII. MULTIPLE PYROPHORIC FLARES 



The main computer program code also simulates the 
dropping of multiple pyrophoric flares, each at different 
instances of time and with generally different parameters. 
These parameters include pyrophoric temperature, T pf , radius 
of pyrophoric at the peak radiant intensity, r^ and the 
time of peak radiant intensity, t . 

Using this option, a scenario can be generated where a 
pre-launched unmanned tactical air launched decoy (TALD) 
drops a sequence of flares in advance of a manned aircraft. 
The missile will then be distracted by a flare generated IR 
noise curtain. The intent is to provide 'cover' to the 
manned aircraft which flies into a safe-zone corridor a few 
seconds later. This is represented in Figure 8.1. 

The settings on the different parameters of the ejected 
pyrophoric can be chosen according to the conclusions 
reached from studying Figures 6.1, 6.2, and 6.4 in Chapter 
VII. In particular, an increase in the pyrophoric 
temperature T pf or an increase in the radius r^ 
significantly decreases S/N between the initial and final 
time of the pyrophoric function. Also, a flare which 
achieves peak radiant intensity at a later point in time 
will have an extended window of effectiveness. 
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The radiant intensity for four pyrophoric flares, which have 
a function burn time of 5s and are launched at times t=0s, 
t=0.5s, t=ls, and t=1.5s are shown in Figure 8.2. The first 
flare has the shortest time to peak radiant intensity, 
t p =0.8s. This value is selected in order to capture the 
missile seeker quickly. Other specifications for the first 
flare include r MAX =2m and T pf =2200K, which were chosen to keep 
S/N low. The second and third flares, which are launched at 
t=0.5s and t=ls, respectively, have a longer time to reach 
peak radiant intensity, t p =1.2s, in order to create an 
extended window of effectiveness. For these two flares, 
r MAX =lm and T pf =2 000K are used. These have been selected lower 
relative to r,^ and T p£ for the first flare since there is 
still significant radiant energy from the first flare. The 
last flare is launched at t=1.5s with larger values of r^ 
and T pf relative to the previous two flares, r MAX =2m and 
T p£ =2200K. Higher values are used in order to maintain a low 
value of S/N by compensating for the reduction in 
effectiveness of the first flare. In addition, t p =1.2s is 
used in order to have an extended window of effectiveness. 

The composite radiant intensity of the four flares 
versus the total function burn time is shown in Figure 8.3. 
The changes in slope are correlated with the activation of 
new flares. 
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The composite S/N versus the total function burn time 
is represented in Figure 8.4. Also, in this plot the changes 
in slope are correlated with the activation of new flares. 

The S/N curves that correspond to each flare are shown 
in Figure 8.5. By comparing Figure 8.4 and Figure 8.5, it is 
obvious there is an extension of the length of time S/N 
remains below some threshold (e.g., 0.5) with the use of 
multiple flares. This point is quantified in Table 8.1. 




ppp p 




Figure 8.1 Pre-launched Unmanned Tactical Air Launched Decoy 
(TALD) drops a Sequence of Flares in advance of a Manned 
Aircraft . 
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J(W/sr) 




Figure 8.2 The Radiant Intensity Distribution of the four 
Pyrophorics . 
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Figure 8.3 The Composite Radiant Intensity Distribution of 
the four Pyrophorics. 



73 



S/N Ratio 




Pyrophoric Function Time (sec) 



Figure 8.4 The Composite S/N versus the Total Function Burn 
Time . 
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Figure 8.5 S/N ratio Curves for Each Flare. 



Table 8.1 Time below tbe 0.5 threshold. 



Flare number 


Time below the 0.5 threshold (sec) 


Flare 1 


~3 


Flare 2 


0 


Flare 3 


0 


Flare 4 


~3 


Composite flare 


~5 . 5 
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In order to demonstrate how these results apply to 
tactical decisions that need to be made in the field, two 
examples are presented. The conditions that the manned 
aircraft be given safe cover while being viewed by a 
"stationary missile" can be stated as 

where 

-At (s/N) is the time S/N is depressed below a specified 
threshold 

-t dae is the approximate time after the dropping of the first 
flare that the aircraft enters the safe corridor 

-V„„ aircraft velocity 

ac 



Equation (VIII-1) can be rearranged as follows: 



V. 



— > 



R< 



^S/ tdac 
/N 



{yin- 2) 



For convenience of interpretation, minimum values for 
(V ac /R f ) are plotted in Figure 8.6 versus At (S/(J) for various 
values in the time delay, t dac . 
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elapsed time S/N depressed 



Figure 8.6 Minimum Velocity Curves for various Values of 
Time Delay Parameter. 



Example 1 

Given V ac =186 m/sec(=0.56 MACH1) , R f =62 m, t dac =1.0s, 
what is the minimum allowed time to keep S/N depressed in 
order to provide "safe" cover? 

Step 1 Calculate V ac /R £ = 3s’ 1 

Step 2 For delay t dac =1.0s, read from Figure 8.6, 1.5s. 

From the above graph and reference to Table 8.1, it is 
obvious that any one of the flares studied would be 
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sufficient. However, if the aircraft delay is changed to 3.0 
seconds and V =62 m/sec, then the condition on minimum At, a ,„, 

changes to about 4.1 seconds and more than one flare packet 
will need to be dropped. 

Example 2 

Given At (s/N) =2s, R £ =62 m, t dac =1.5s, 
what is the minimum velocity of the aircraft in order to be 
provided "safe" cover? 

Step 1 For At (S/N) =2s and t dac =1.5s, read from Figure 8.6, 

(V ac /R £ ) mi ~ 2.0s' 1 

Step 2 Use the results from step 1 with R £ =62 m to obtain 
(V ac)MIN =124 m/sec (=37% MACHl ) 

The point of this example is to demonstrate that a time 
limit in the required reduction of S/N forces a predictable 
lower limit on the velocity of the manned aircraft. 

From Figures 8.4 and 8.6 several significant and 
general observations can be deduced. First, because the 
period of time in which S/N stays below a given threshold is 
finite, as can been seen from Figure 8.4, there is an upper 
limit on the delay that can take place before the manned 
aircraft enters the safe corridor. The longest this delay 
can be is the total time interval that S/N is below the 
threshold. Note from Figure 8.6 that each curve predicts the 
minimum required velocity of aircraft which approaches 
infinity as the elapsed time for S/N to be depressed 
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asymptotically approaches the delay. For example, the curve 
with a delay of three seconds approaches infinity from the 
right as the elapsed time S/N depressed approaches three 
seconds . 

These curves predict in general that a greater delay in 
the aircraft entering the corridor requires a greater 
velocity of the manned aircraft in order to guarantee safe 
passage. Conversely, the curves when applied in reverse 
suggest a minimum required time for the depression in S/N. 
In particular, from Figure 8.6 it follows that the greater 
the delay and/or the slower the velocity of the manned 
aircraft, the longer S/N needs to stay depressed below a 
specified threshold. 



79 



80 



IX. CONCLUSIONS AND FUTURE ENHANCEMENTS 



The objective of this research was to model the impact 
of dropping pyrophoric type flares which are in the field- 
of-view of a prelaunch missile. It is assumed that these 
flares are launched from an unmanned tactical air launched 
decoy (TALD) . The manned aircraft, which needs protective 
cover and will presumably benefit from the flares, will 
follow into the corridor a few seconds after the start of 
the flare sequence. 

In order to reach this objective, a simulation model of 
the time-space radiant intensity distribution for a 
pyrophoric flare expendable has been developed. The 
performance of the pyrophoric flares to create a distraction 
for the missile was characterized in term of S/N where S is 
the radiant power of the plume and N is the radiant power of 
one or more flares. Several conclusions from this study are 
described in the next two paragraphs . 

The radiant power for the pyrophoric flare burn in the 
[3-5] micron band was found to be significantly higher than 
what occurs in the [8-12] micron band. It was also found 
that increasing the pyrophoric temperature T £ or increasing 
the radius r,^ significantly decreases S/N between the 
initial and final time of the pyrophoric function. Finally, 
a flare which achieves a peak radiant intensity later in 
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time has an extended window of effectiveness; i.e., a longer 
window in time for which S/N remains below a specified 
threshold . 

In this thesis the concept of dropping of multiple 
pyrophoric flares to create a "safe" corridor for manned 
aircraft missions was evaluated. In particular, multiple 
flare "packets" lengthen the window over which S/N remains 
below an acceptable threshold. 

Lastly, it was noted that if the manned aircraft is not 
moving with high enough velocity there will be a problem in 
generating "safe" cover. From this concept, a rule relating 
minimum aircraft velocity to both the delay in entering the 
safe corridor and the time S/N is depressed was presented 
and illustrated with several examples. 

A number of simplifying assumptions were applied in the 
study. Future work on the simulation model could focus on 
one or more of these assumptions. For example, to mention a 
few points, atmospheric absorption was neglected in this 
study. Also, the plume was assumed to be a constant 
temperature graybody. A better model for the atmospheric 
absorption would be to consider the cold C0 2 absorption 
spectrum. For the plume, consideration of the hot C0 2 
emission spectrum would be also beneficial (the red spike- 
blue spike effect) . Also, an experimental data based model 
for the pyrophoric burn radiation spectrum could be used to 
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refine the simulation model. Finally, the calculations of 
S/N could be based on more complex image processing 
algorithms that are indicative of modern reticle based 
missiles . 
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APPENDIX A. 



COMPUTER AND ANALYSIS VARIABLES 



Computer 

Variables 


Analysis 

Variables 


Brief Description 


FOV 


FOV 


Missile f ield-of-view 


R 


R 


Range between missile and target 


Rf ield 


^Field 


Range inside missile's f ield-of- 
view on target's plane 


Tpf 


T 

1 pf 


Pyrophoric temperature 


Tpl 




Plume temperature 


lambda 1 


*x 


Initial wavelength 


lambda 2 




Final wavelength 


tl 


tx 


Initial time of pyrophoric 
function 


t2 


t 2 


Final time of pyrophoric function 


t_peak 


fc p 


Pyrophoric time at max radiant 
intensity 


r_max_t 


T 

MAXt 


Radius of pyrophoric 
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rmax 


2T 

MAX 


Radius of pyrophoric at max 
radiant intensity 


rnorm 


NORM 


Normalization radius of 
pyrophoric 


Afl 


A fl 


Area of pyrophoric at max radiant 
intensity 


Apl 


A pi 


Area of plume 


a 


a pi 


Major axis of plume's ellipse 
area 


b 


b pi 


Minor axis of plume's ellipse 
area 


vel 


U p 


Velocity of pyrophoric 


e 


e 


Average emissivity 


s 


o 


Stef an-Boltzmann constant 


Wf lare 


W(X,T) pf 


Total emittance of pyrophoric 


Wplume 


W(X,T) pl 


Total emittance of plume 


Wpfb 


w(^A 2 ) p£ 


Emittance of pyrophoric in the 
band of interest 


Wplb 


w(X,X 2 ) pl 


Emittance of plume in the band of 
interest 


ratWf 1 




Fraction of pyrophoric emittance 
in the band of interest 
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ratWpl 


*0 

H* 


Fraction of plume emittance in 
the band of interest 


J_max 


^MAX 


Peak radiant intensity of 
pyrophoric 


J 


J(t) 


Radiant intensity of pyrophoric 
function 


Jpfb 




Radiant intensity of pyrophoric 
in the band of interest 


Jplb 




Radiant intensity of plume in the 
band of interest 


fa 


f (t,a) 


Gamma distribution 


G 


r (a) 


Gamma function 


a 


a 


Parameter (a) of Gamma 
Distribution 


fad 


f G) 


Gaussian space distribution 


s 


a 


Standard deviation of Gaussian 
distribution 


m 


n 


Mean of Gaussian distribution 


Q 


Q 


Amount under a Gaussian 
distribution curve 


scrJ 


j (JW t/t p ) 


Weighted value of a Gaussian 
distribution in normalized space 
domain 


Intscr 


J v (r 0 , t) 


Integrated view radiant intensity 


Npix 


N 'pix 


Number of pixels 


Cpf 


c« 


Pyrophoric image matrix 
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Cpl 


c « 


Plume image matrix 


Cpfn 


C'„ 


Normalized pyrophoric image 
matrix 


Cpln 


c 'ii 


Normalized plume image matrix 


Ratio 


S/N 


S/N ratio 
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APPENDIX B. MATLAB CODES 

This appendix contains listing of all MATLAB input 
files that were used to get the results posted in Chapters 
VI, VII and VIII. These M-files are: 

1. Pyrof.m, 

2 . Plankpf .m 

3 . Plankpl .m 

4 . Dquadl .m 

5. Dquadl. m, and 

6 . Gamma z . m 

Pyrof.m is the main computer code which has two parts. 
The first part generates the S/N curves for one or multiple 
pyrophoric with different parameters and the second part 
generates the pyrophoric and plume images . 



Plankpf. m is a function file which calculates the 
Spectral Radiant Emittance (SRE) for the pyrophoric at 
temperature T pf (Kelvin) . 



Plankpl. m is 


a function 


file 


which 


calculates 


the 


Spectral Radiant 


Emittance 


(SRE) 


for 


the 


plume 


at 


temperature T pl (Kelvin) . 












Dquadl. m is 


a function 


file 


which 


calculates 


the 


Spectral Radiant 


Emittance 


(SRE) , 


for 


the 


plume 


at 



temperature T pl (Kelvin), inside the detectors band. 



89 



Dquad2.m is a function file which calculates the 
Spectral Radiant Emittance (SRE) , for the pyrophoric at 
temperature T p£ (Kelvin), inside the detectors band. 

Gammaz.m is a function file which calculates the GAMMA, 
function . 
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Program 1 Pyrof . m 



%%%%%%%%%% MAIN COMPUTER CODE %%%%%%%%%% 
%%%%Ref. files: 

% plankpf.m, plankpl.m, 

% dquadl.m, dquad2.m, gamxnaz.m 



format long 

global Tpl 
global lambdal 
global lambda2 
global Tpf 
global tl 
global t2 

%%%%% Inputs %%%% 

program=input ( ' For the S/N part enter 1, else for the images part enter 
2 : ' ) ; 

%%%%%%%%%%%%%%%%%%%%% First Part 
if program==l 

Tpl=input ( ' Set the Temperature of plume (Kelvin): '); 

Tp=input ( ' Set the Temperature of Pyrophoric (Kelvin): '); 

lambdal= input { ' Set the value of the lambdal (um) : ' ); 
lambda2=input ( ' Set the value of the lambda2 (um) : ' ) ; 

tl=input { ' Set the initial value of the time for the Pyrophoric (sec): 

') ; 

t2=input ( ' Set the final value of the time for the Pyrophoric (sec ) : '); 

t_peakl=input ( ' Set the value of the t_peak (sec): '); 

rmaxp=input ( ' Set the value of the rmax (m) : '); 

t3=input ( ' Set the increment step between tl and t2 : '); 

tf l=input ( ' Set the dropped times of pyrof oric : '); 

t=tl : t3 : t2 ; 
figure (1) 

for z=l:length(tfl); 
tf lare=tf 1 ( z ) ; 
if length (t_peakl ) ==1 ; 

t_peak=t_peakl ; 
else 

t_peak=t_peakl ( z ) ; 

end 

if length (Tp) ==1 ; 

Tpf =Tp; 
else 

Tpf =Tp ( z ) ; 

end 

if length (rmaxp) ==1 ; 

rmax=rmaxp; 

else 

rmax=rmaxp ( z ) ; 

end 

e=0.9; %emissivity of Pyrophoric 

s=5 . 67*10 A ( -12 ) ; % Stef an-Boltzmann constant (cm) 
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Wpf b=dquad2 ( ' plankpf 7 , lambdal , lambda2 ) ; %Emittance between the two lambda 
%( W cm A ( -2 ) ) 

Wf lare=s*Tpf A 4 ; %Stefan-Boltzmann ( total emittance over all wavelengths) 
Af l=4*pi* ( rmax) A 2 ; % Area Pyrophoric at t_peak 
sl = 5 . 67*10 A ( -8 ) ; % Stef an-Bol t zmann constant (m) 

J_max= ( (e*sl*Tpf A 4) / (4*pi) ) *Af 1 ; % Total Radiant Intensity (W/sr) 
ratWf 1= (Wpfb/Wf lare) ; 

m= ( tf lare/ 1 3 ) +1; 
a=t_peak+l ; 

G=gammaz (a ) ; % Find 'G' with the help of 'gammaz 7 function 
fa= ( t . A ( t_peak) . *exp ( -t ) ) . /G; % gamma equation 

J= ( f a . * J_max) * ( 1/max ( fa) ) ; % Normalized to the total Radiant Intensity 
% J_max 
Js=J*ratWf 1 ; 

[rl , r2 ] =size ( Js) ; 
cl=Js ( 1 , 2 : r2 ) ; 
c= [zeros (l,m) ,cl] ; 
if z = = 1 ; 

y=l : length (c) ; 
y=c; 

else 

y= [y , zeros ( 1, ( length (c ) -length (y) ) ) ] +c; 

end 

tpl=0 : t3 : ( ( length (c) -1) *t3 ) ; 

eval ( [ 7 save data 7 , num2str (z ) , 7 . mat ' ] ) ; 

eval ( [ 7 load data 7 , num2str (z) , ' .mat 7 ]); 

plot ( tpl / c ) 

grid on 

xlabel ( 7 time (sec) 7 ) 
ylabel ( 7 J (W/sr) 7 ) 
hold on 
end 

figure (2) 

tcomp=0 : t3 : ( ( length (y) -1) *t3 ) ; 
plot ( tcomp, y) 
grid on 

xlabel ( 7 time (sec) 7 ) 
ylabel ( 7 J (W/sr) 7 ) 

% Calculation of radiant intensity for the plume 

Wplb=dquadl ( 7 plankpl 7 , lambdal, lambda2) ;%Emittance between the two lambda 

%( W cm A ( -2 ) ) 

a=1.50; %minor axis (m) 

b=3.50; %major axis (m) 

Apl=pi*a*b; % (Area) 
e=0 . 9 ; %emissivity 

s=5 . 67*10 A ( -12) ; %Stef an-Bol tzmann constant 

Wplume=s*Tpl A 4 ; %Stefan-Bolt zmann ( total emittance over all wavelengths) 

sl=5 . 67*10 A ( -8 ) ; % Stef an-Bol tzmann constant (m) 

ratWpl=Wplb/Wplume; 

Jpl= ( (e*sl*Tpl A 4) /pi) *Apl*ratWpl; %Radiant Intensity between the two 
%lambda (W/sr) 

Jplume=Jpl . *ones ( 1 , length (y ) ) ; 

^ ********* 5 /i^ tio ******* 

Ratio=Jplume . /y ; 

figure (3 ) 

plot ( tcomp , Ratio) 
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grid on 

xlabel ( ' Pyrophoric Function Time (sec) ') 
ylabel('S/N Ratio') 

figure (4) 

for x=l : length ( tf 1) ; 
tf lare=tf 1 (x) ; 
if length ( t_peakl) = = 1 ; 

t_peak== t_peakl ; 
else 

t_peak=t__peakl (x) ; 

end 

if length (Tp) ==1 ; 

Tpf ==Tp ; 
else 

Tpf =Tp (x) ; 

end 

if length (rmaxp) ==1 ; 

rmax==rmaxp; 

else 

rmax=rmaxp (x) ; 

end 

% Calculation of radiant intensity for the pyrophoric 

e=0.9; %emissivity of Pyrophoric 

s=5 . 67*10 A (-12) ; % Stefan-Boltzmann constant (cm) 

Wpfb=dquad2 ( 'plankpf ' , lambdal, lambda2) ; %Emittance between the two lambda 
%( W cm A ( -2 ) ) 

Wf lare=s*Tpf A 4 ; %Stefan-Boltzmann ( total emittance over all wavelengths) 
Afl=4*pi* (rmax) A 2; % Area Pyrophoric at t_peak 
sl=5 . 67*10 A (-8) ; % Stefan-Boltzmann constant (m) 

J_max= ( (e*sl*Tpf A 4) / (4*pi) ) *Af 1; % Total Radiant Intensity (W/sr) 
ratWf 1= (Wpfb/Wf lare) ; 

m= ( tf lare/ 13 ) +1 ; 
a=t_peak+l ; 

G=gammaz (a) ; % Find 'G' with the help of 'gammaz' function 
fa= (t . A (t_peak) . *exp(-t) ) . /G; % gamma equation 

J=(fa. *J_max) * (1/max (fa) ) ; % Normalized to the total Radiant Intensity 
%J_max 

Js=J*ratWfl; 

[rl, r2] =size ( Js) ; 
cl=Js ( 1 , 2 : r2 ) ; 
c2= [zeros (l,m) , cl] ; 

Jplume=Jpl . *ones ( 1 , length (c2 ) ) ; 
c=Jplume . / c2 ; 

tpl=0 : 1 3 : ( ( length (c) -1) *t3) ; 

eval ( [ ' save datas ' , num2str (x) , ' .mat ' ] ) ; 
eval ( [ ' load datas ' , num2str (x) , ' .mat ' ] ) ; 
plot ( tpl, c) 
grid on 

xlabel ('time (sec)') 
ylabel('S/N Ratio') 
hold on 
end 

%%%%%%%%%%%%%%%%%%%%% Second Part 
else 

Npix=input ( ' Set the number of pixels (it must be an integer value with 
half an odd number) : ' ) ; 

tspf =input ( ' Set the time (sec) that you want to create the images: '); 
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Tpl=input ( ' Set the Temperature of plume (Kelvin): '); 

Tpf =input ( ' Set the Temperature of Pyrophoric (Kelvin): '); 

lambdal= input (' Set the value of the lambdal (um) : '); 
lambda2=input ( ' Set the value of the lambda2 (um) : '); 

tl=input ( ' Set the initial value of the time for the Pyrophoric (sec): 
'); 

t2=input ( ' Set the final value of the time for the Pyrophoric ( sec ) : '); 
t__peak=input ( ' Set the value of the t__peak (sec): '); 

rmax=input ( ' Set the value of the rmax (m) : '); 

t3=input ( ' Set the increment step between tl and t2 : '); 

% Calculation of radiant intensity for the plume 
Wp.lb=dquadl ( 'plankpl ' , lambdal, lambda2) ; 
a=1.5; %minor axis (m) 
b=3.5; %major axis (m) 

Apl=pi *a*b; % (Area) 
e=0 . 9 ; %emissivity 

s=5 . 67*10^ ( -12 ) ; %Stef an-Boltzmann constant 
Wplume=s*Tpl A 4 ; %total emittance over all wavelengths 
sl=5 . 67 *10 A ( -8 ) ; % Stef an-Boltzmann constant (m) 

Jplb= ( (e*sl*Tpl / "4) /pi) *Apl* (Wplb/Wplume ) ; %Radiant Intensity between the 
%two lambda (W/sr) 

F0V= ( (pi*2) /180) ; 

Range=l . 8* ( 10 A 3 ) ; % (m) 

Rfield=FOV* Range; % (m) 

%I create the image matrix for the Plume 
Cpl=zeros (Npix+1 , Npix+1) ; 
for il=-Npix/2 : 1 :Npix/2 ; 
for jl=-Npix/2 : 1 : Npix/2 ; 
ipixl=il+ ( (Npix/2 ) + 1) ; 
jpixl= j 1+ ( (Npix/2) +1) ; 

rl= ( (il) . ^2 . / (a/Rf ield) ^2 ) + ( (jl) . A 2 . / (b/Rf ield) A 2 ) ; 
if rl<=(Npix+l)^2; 

Cpl ( ipixl, jpixl) =Jplb; % I assign the value Jplb to the pixels 
%inside the matrix 
else 

Cpl (ipixl, jpixl) =0;% I assign the value 0 to the pixels outside 
%the matrix 

end 

end 

end 

suml = 0 ; 

for i=l : (Npix+1 ) ; 
for j =1 : (Npix+1 ) ; 

suml= suml+ Cpl(i,j); 

end 

end 

newCpl= Cpl/suml *Jplb; 

Cpln=newCpl ; 

% Calculation of radiant intensity for the pyrophoric 

e=0.9; %emissivity of Pyrophoric 

s=5 . 67*10^ (-12) ; % Stef an-Boltzmann constant (cm) 

Wpfb=dquad2 ( 'plankpf ', lambdal, lambda2 ); %Emittance between the two lambda 
%( W cm A (-2) ) 

Wf lare=s*Tpf ^4 ; %total emittance over all wavelengths 
Af 1=4 *pi* (rmax) A 2 ; % Area Pyrophoric at t_peak 
sl=5 . 67*10^ (-8) ;% Stef an-Boltzmann constant (m) 
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J_max= ( (e*sl*Tpf "M ) / ( 4*pi) ) *Af 1 ; % Total Radiant Intensity (W/sr) 
ratWf 1= (Wpfb/Wf lare) ; 

t=tl : 0 . 1 : t2 ; 
vel=rmax. /t_peak; 
r_max_t=t . *vel ; 
rnorm= r_max_t / rmax ; 

% I create the time distribution (gamma distribution) for a=t_peak+l 
a=t_peak+l; 

G=gammaz (a) ; % Find ' G ' with the help of ' gamma z ' function 
f a= ( t . * ( t_peak) . *exp ( -t ) ) . /G; % gamma equation 

J= ( f a . * J_max) * ( 1/max ( fa) ) ; % Normalized to the total Radiant Intensity 
% J_max 



% I create the normalized space distribution 
s=l ; % standard deviation 

for tx=tl : t3 : t2 ; 

f ad= (exp ( - (rnorm- ( tx/t_peak) ) . A 2 . / (2*s A 2) ) ./ (sqrt (2 *pi) *s) ) ; %Gaussian 
%distribution 
12=tx; 

f a2= ( 12 A ( t_peak) *exp ( -12 ) ) /gammaz ( t_peak+l ) ; 

J2= (fa2*J_max) * (l/max(fa) ) ; % Normalized to the total Radiant 
% Intensity J_max 

Q2= ( 0 . 5) . *erfc ( ( -tx/ (t_peak*s) ) . /sqrt (2) ) ; % Integral from zero to 
%infinity for a Gaussian distribution with mean at tx/t_peak 
scrJl=f ad*J2/Q2 ; %The weighted values for every Gaussian distribution 
%with mean at tx/t_peak 
end 

for z=l : length ( tspf) ; 
ts=tspf ( z) / 

fay= ( exp ( - ( rnorm- ( ( ts/t_peak) )) .~2./(2*s /N 2)) ./( sqrt ( 2 *pi) *s) ) ; 
ls=ts ; 

fas= ( ls A ( t_peak) *exp ( -Is) ) /gammaz ( t_peak+l ) ; 

Js= ( fas* J_max) * (1/max ( fa) ) ; 

Qs= (0.5) . *erf c ( ( -ts/ ( t_peak*s) ) . /sqrt (2 ) ) ; 
scrJ3=fay* Js/Qs; 

%I create the integrated Pyrophoric model 
Np=3 0 ; 

[rl # cl ] =size (rnorm) ; 
romax=max (rnorm) ; 

Nro=cl-l ; 

Intscr2=l :Nro+l; 

ro=l:Nro+l; % ro is the distance from the center of projection 
Dro=romax/Nro ; 
for J=1 :Nro+l ; 
sum=0; 

ro ( J) = ( J-l ) . *Dro; % for every J pick up one ro 
xmax=sqrt ( (romax) . A 2- (ro ( J) ) .* 2 ) ; 

Dx=xmax. /Np; 
for n=l:Np+l; 
x= (n-1 ) . *Dx ; 

r=sqrt (x.^2+(ro(J) ) . A 2) ; 
kl= (r . /max (rnorm) ) *cl; 

k=round(kl); % index k goes to closest integer 
if k>cl 
a=0; 

sum=sum+a; 
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elseif k==0; 
a=0 ; 

sum=sum+a ; 
else 

sum=sum+scrJ3 (k) ; 
end 

end 

Intscr2 ( J) =sum; 

end 

%I create the image matrix for the Pyrophoric model 
FOV= ( (pi*2)/180) ; 

Range=l . 8* (10^3 ) ; % (m) 

Rf ield=FOV*Range; % (m) 

[rl , cl ] =size ( rnorm) ; % cl takes the value of size ( r_max_t ) 
Cpf=zeros (Npix+1 , Npix+1 ) ; 
for i=-Npix/2 : 1 :Npix/2 ; 
for j = -Npix/2 : 1 :Npix/2 ; 
ipix=i+ ( (Npix/2 ) +1 ) ; 
jpix= j + ( (Npix/2) +1) ; 

Nr=sqrt ( i . ^2+ j . A 2 ) ; 

R= (Rf ield/ (Npix+1 ) ) *Nr ; 

Ll= (R/ (max ( rnorm) *rmax) ) *cl ; 

L=round(Ll) ; 
if L>cl 

Cpf (ipix, jpix) =0; 
elseif L-=0 ; 

Cpf ( ipix, jpix) =0; 
else 

Cpf (ipix, jpix) =Intscr2 (L) ; 

end 

end 

end 

Jts=Js*ratWf 1; 
sum=0 ; 

for i=l : (Npix+1 ) ; 
for j=l : (Npix+1) ; 

sum= sum+ Cpf(i,j); 

end 

end 

newCpf= Cpf/sum *Jts; 

Cpfn=newCpf ; 



% ********* s/N Ratio ******* 

Ratio=Jplb/ Jts 

eval ( [ # save data ' # num2str (z) , ' .mat ' ] ) ; 
end 

for z=l : length (tspf) ; 

eval ( [ x load data' , num2str {z) , ' .mat ' ] ) ; 
figure (z) 

image (Cpfn, r CDataMapping y # 'scaled') 
end 

figure ( z+1 ) 

image (Cpln, ' CDataMapping ' , 'scaled') 
end 
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Program 2 Plankpf.m 



% This function file calculates the Spectral Radiant Emittance (SRE) 
% for the pyrophoric at temperature Tpf (Kelvin) . 

function Y = plankpf (Lambda, T) 
global Tpf 
cl=3 .7415e4; 
c2=l . 4388e4 ; 

Y = cl ./(Lambda .^5 . * (exp(c2 ./(Lambda . * Tpf))-1)); 

% Y will have units of watts/ (cm^2-micron) 



Program 3 Plankpl.m 

% This function file calculates the Spectral Radiant Emittance (SRE) 
% for the plume at temperature Tpl (Kelvin). 

function Y = plankpl ( Lambda , T ) 
global Tpl 
cl=3 .7415e4; 
c2 = l . 43 88e4 ; 

Y = cl ./(Lambda . A 5 . * (exp(c2 ./(Lambda . * Tpl))-1)); 

% Y will have units of watts/ (cm^2-micron) 



Program 4 Dquadl . m 

% DQUADl Numerically integrates an expression using QUAD 8 . 

function [I] =dquadl (Fun,al, a2 ,a3 ,a4) 

var = ' x ' ; 
tol=le-3 ; 
if nargin<3, 
help dquadl 
return 
end 

if nargin == 3, 
xmin=al; 
xmax=a2 ; 

elseif nargin == 4, 
if (isstr (al) ) , 
var=al ; 
xmin=a2; 
xmax=a3 ; 
else 
xmin=al ; 
xmax=a2 ; 
tol=a3 ; 
end 
else 

var=al ; 
xmin=a2 ; 
xmax=a3 ; 
tol=a4 ; 

end 

I=quad8 ( 'plankpl ' / xmin / xmax / tol) ; 
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Program 5 Dquad2 . m 



% DQUAD2 Numerically integrates an expression using QUAD 8 . 

function [I] =dguad2 (Fun, al , a2 , a3 , a4) 

var = 'x ' ; 
tol=le-3 ; 
if nargin<3, 
help dquad2 
return 

end 

if nargin == 3, 
xmin=al ; 
xmax=a2 ; 

elseif nargin == 4, 
if (isstr (al) ) , 
var=al ; 
xmin=a2 ; 
xmax=a3 ; 
else 

xmin=al ; 
xmax=a2 ; 
tol=a3 ; 
end 
else 

var=al; 
xmin=a2 ; 
xmax=a3 ; 
tol=a4 ; 

end 

I=quad8 ( 'plankpf ' , xmin, xmax, tol) ; 



Program 6 Gamma z .m 

% GAMMAZ gamma function that also allows for complex arguments, unlike 
% MATLAB's GAMMA function. For real arguments, GAMMAZ gives the 

% same results as GAMMA (within numerical error) . 

function [f] = gammaz(z); 

[m n] = size ( z ) ; 
f = zeros (m, n) ; 
cof = [76.18009172947146 
-86.50532032941677 
24.01409824083091 
-1.231739572450155 
0.1208650973866179e-2 

-0 . 5395239384953e-5] ; % 6 coefficients in series expansion 

for icol = l:n, % do one column at a time 

zz = z ( : , icol) ; 

zp = ones ( 6 , 1 ) *zz ' + [ 1 : length (cof )]' *ones ( 1 , length ( zz )) ; % vectorize 
ser = (cof '* (l./zp)+l. 000000000190015) ' ; 
tmp - zz + 5.5 - ( zz+ . 5) . *log (zz + 5 . 5 ) ; 

Ingamma = -tmp + log(2.5066282746310005*ser./zz); 
f ( : , icol ) = exp ( Ingamma) ; 
end 
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APPENDIX C. 



MATHEMATICAL DERIVATIONS 



1. Solution for t =a-l 

P 

The pdf of a standard Gamma RV is given by equation (III- 
7) in Chapter III. For convenience it is reproduced below 

(c-i) 

t > 0 
a > 0 



Now by setting 

4:fM - o 

at 

with r(a) = (a - 1)!^ 0 for a > 0 



we get 



rW 




+ l* _l (-l)e'']=0=> 



==> [(£r — l)r a 2 e 1 +t a ' - {— l)e 'J=0=> a = t + l 



(C-2) 



(C-3) 



Now, for t=t 



t 



p 



a - 1 



(C-4) 
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2 . Positive side integration of Gaussian distribution 

The pdf of a Gaussian RV is given by equation (III-10) in 
Chapter III. For convenience it is reproduced below 



r) 



1 

V2^r -cr 



Jx-M) 2 / 

e 



— oo < X < °° 

- oo < fj. < oo mean 

0 < <7 5 tan dard deviation 



(C-5) 



The above equation for |A=0 and a=l becomes 







— oo < X < oo 



(C-6) 



To find the area under the curve from x to infinity the 
following integral must be evaluated: 



G(*;0,1) = 



l 

V2-7T 



r 




(C-7) 



In terms of the complementary error function 



Q{x) = Y 2 -erf c 

Now, if the continuous random variable X has mean |J. and 
standard deviation o then 
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(C-9) 



z = 



X-LI 



a 



is a Gaussian random variable with ^.=0 and G=l, and the area 
under the curve from z to infinity is 



Q(z; 0,1) 




(C-10) 



From equation (C-8), we have 




(C-ll) 



The area under the curve from zero to infinity is obtained 
by setting x=0 in the previous equation to obtain [5] 






(C- 12) 
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